clip-02-05

clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
gr(size=(600,300));
Plots.GRBackend()

snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]
5-element Array{Float64,1}:
 0.881
 0.64 
 0.818
 0.632
 0.778

snippet 3.3

Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);
10000-element Array{Float64,1}:
 0.493
 0.839
 0.571
 0.919
 0.7  
 0.815
 0.589
 0.549
 0.578
 0.604
 ⋮    
 0.572
 0.543
 0.701
 0.535
 0.686
 0.779
 0.628
 0.716
 0.814

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);
Object of type Chains, with data of type 10000×1×1 Array{Float64,3}

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

parameters
      Mean    SD   Naive SE  MCSE     ESS  
toss 0.6357 0.1407   0.0014 0.0015 8285.226

Describe the chain

describe(chn)
Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

Empirical Posterior Estimates
───────────────────────────────────────────
parameters
      Mean    SD   Naive SE  MCSE     ESS
toss 0.6357 0.1407   0.0014 0.0015 8285.226

Quantiles
───────────────────────────────────────────
parameters
      2.5% 25.0% 50.0% 75.0% 97.5%
toss 0.061  0.54 0.644 0.741 0.984

Plot the chain

plot(chn)
0 2500 5000 7500 10000 0.2 0.4 0.6 0.8 1.0 toss Iteration Sample value 0.00 0.25 0.50 0.75 1.00 0.0 0.5 1.0 1.5 2.0 2.5 toss Sample value Density

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws

snippet 3.5

Analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

Add quadratic approximation

plot(p..., layout=(1, 2))
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws 0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

End of 03/clip-02-05.jl

This page was generated using Literate.jl.